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Abstract 

Recently semiclassical approximations have been successfully applied 
to study the effect of a superconducting lead on the density of states and 
conductance in ballistic billiards. However, the summation over classical 
trajectories involved in such theories was carried out using the intuitive 
picture of Andreev reflection rather than the semiclassical reasoning. We 
propose a method to calculate the semiclassical sums which allows us to go 
beyond the diagonal approximation in these problems. In particular, we 
address the question of whether the off-diagonal corrections could explain 
the gap in the density of states of a chaotic Andreev billiard. 

PACS; 05.45.Mt, 74.45.+C 

1 Introduction 

Andreev billiards are two-dimensional electron cavities with a superconduct- 
ing part of the boundary (see [1] for a review) . A negatively charged electron 
incident on the superconductor is reflected as a positively charged hole with 
an opposite momentum. This process of Andreev reflection [2] generates a 
new kind of dynamics in comparison to conventional (normal) billiards. The 
excitation spectrum of an Andreev billiard depends on the shape of its bound- 
ary [3,4]. Using the random-matrix theory, it was shown [3] that the density 
of states (DoS) d{E) in the chaotic billiard has a minigap around the Fermi 
energy Ep-, while in the integrable billiard it is proportional to E, the energy 
counted from Ep- The width of the minigap is of the order of the Thouless 
energy, which is much smaller than the superconducting gap A. The semiclas- 
sical theory of [4] confirms the above results in the integrable case, but finds an 
exponential suppression of d{E), instead of a gap, in the chaotic cavity. The 
authors surmised that the disagreement can be attributed to the use of diag- 
onal approximation. The method proposed below allows one to determine the 
off-diagonal corrections in a billiard with a superconducting lead. (A somewhat 
similar concept was used in [5] in the case where the whole boundary is super- 
conducting.) We find that these corrections may, indeed, reflect the existence 
of the gap, although the full understanding of the problem is still missing. 
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Another area of the apphcation of the semiclassical techniques is transport 
problems. In general, the conductance of a billiard with a superconducting 
lead can be expressed in terms of transmission and reflection amplitudes with 
or without clcctron-to-hole conversion [6,7]. In particular, the contribution of 
the Andrecv reflection to the conductance is proportional to Ra = Tr (rher^^), 
where rhe is a part of the scattering matrix describing an electron injected into 
the billiard returning back as a hole. In [8], this quantity was calculated semi- 
classically under the assumption of the exact Andreev reflection {E — > 0) for a 
three-probe geometry. Both the diagonal and off-diagonal (weak localisation) 
contributions were given. It was argued there that Ra is equal simply to the 
normal transmission coefficient for an electron to reach the superconducting 
lead (then, classically, it would return back as a hole with a 100% probabil- 
ity). Hence, the off-diagonal corrections to Ra could be computed by standard 
methods [9] developed for a normal billiard. In this letter, we start from the 
quantum relation [10] between Ra and the normal transmission amplitudes (in 
the case of a two-probe geometry) and then apply the semiclassical limit. Thus, 
no additional classical assumption is required. It is found that, while the di- 
agonal part of Ra is, in fact, equal to the classical normal transmission, the 
weak-localisation correction is nontrivial. 



2 Density of states 

2.1 Problems with the ecirUer approach 

Here, we outline the semiclassical theory of [4] and point out the technical 

difficulties that arise when sums over classical trajectories arc calculated. We 
consider a chaotic billiard with a superconducting lead. The lead width w is 
assumed to be small compared to the perimeter of the billiard, so as not to 
spoil its chaoticity. The energy is taken to be high enough to provide for a large 
number of open channels N = kFw/ir ^ 1 in the lead {kp is the Fermi wave 
number). The lead is modelled as an ideal wire connecting the billiard and the 
superconductor. The average quasiparticle DoS is given by [4] 

d{E) = do - -Im^ ^Y^— (Tr [S{E)S* {-E)]') , (1) 

where do = Mk/2-Kfi? is the average DoS for a particle of mass M in the isolated 
billiard of area A, S{E) is the electron scattering matrix for the billiard with 
the normal lead, the energy E is measured from Ep and the average is taken 
over a classically small interval oikp- In the semiclassical approximation, the 
scattering matrix [11, 12] 

Snm{E)= J2 A^expUsA (2) 

7(n,m) 

is a sum over classical trajectories j{n,m) starting and ending at the lead 
and connecting the channels m and n. The trajectories make an angle 9„i = 
± sm~^ {rmr/kFw) with the lead direction as they enter the billiard, and 9n — as 
they leave it. Here we assumed \E\ <C Ep. The actions then can be expanded 
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about Ep as Sj{E) ~ «S-y(0) + ETy{0), where Ty is the time and 

S.y{pn,Pm;E = 0) = / P ■ dr - PnVn + PrnVm 



(3) 



is in the momentum representation on the lead. Here, ym{n) are initial (final) 
coordinates at the lead cross section and Pm{n) = fikp sin 9m{n) are the conjugate 
momenta. The prefactors A-y = \A-y\ exp(— i/U-y7r/2) depend on the the Maslov 
indices jj^ and 



2kFW^ 



dVn 



9(sin ( 



1/2 



2^ 



dpndp„ 



1/2 



(4) 



Using the semiclassical expression (2), the authors of [4] arrive, in the diag- 
onal approximation, at the following result: 



N 



{T:T[S{E)S*{-E)f)= ^ J2 l^l'expU2/£;TJ 

n,m=l 7(n,m) 



2^ 



d{sm9) / dyexp 
1 Jo 



h 



2lET{y, 



(5) 



where the sums become the integrals over the initial (or final) conditions (y, 6) 
in the limit N ^ 1. While the derivation of this formula for i = 1 is given in [4] , 
its generalization to Z > 1 is not obvious. For example, when I = 2,we can write 

TV [5(^)5* (-£)]2 



JV 



= E 



E 



•^71 '^7( •^72 -^7^ exp 



m,"!, 7i(ni,n'i),7j(n'i,n2), 

"2, "2 = 1 72(n2,ra2)>72("2>"l) 



i^Ji ~ '^7i + ^^72 ~ '^72) 



X exp 



-E {Xf^ +Ty^+T^^+ Ty^) 



(6) 



After this expression is averaged over a small window of kp, only the terms where 
the actions cancel in the exponent will remain. The diagonal approximation 
amounts to keeping the terms with all four actions equal. In a chaotic system 
without symmetries this is only possible when 71 = {'y'i)~^ = 72 = (72)^^; 
where 7"^ is the time reversed path 7. This condition corresponds to the exact 
Andreev reflection at £^ = 0. It yields the average 



{Tr[S{E)S*{-E)]') = J2 E l-47l'exp(-4£;T7 

n,n' 7(71.7?,') 



(7) 



which is a wrong result. Thus, this calculation shows that the diagonality 
requirement alone is not sufficient in the case of Andreev billiards. As we shall 
see next, the stationary-phase (SP) approximation must be employed when 
calculating the semiclassical sums. 



2.2 Traces of scattering matrices 

In this section, we derive equation (5) and find the quantum corrections to it. 

In order to focus on the idea of the method, the case Z = 2 is considered here 
and the general situation is treated in appendix A.l. 
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2.2.1 Classical contribution 

We begin by rewriting tlie cliannel sum (6) as an integral over momenta 



T:v[S{E)S*{-E)f =( — ) / dpidp[dp2dp'. 



'■'2 X! 

72,72 



A A*, A A*, 

•^71 •^7( -^72 -^7^ 



X exp 



exp 



-E {Ty^ + TY^+Ty^+ Ty^ 



(8) 



Although each channel allows for two signs of the momentum, it is expected 
that the leading contribution to the integral comes from the Andreev-reflected 

paths. Therefore, given pi, p'l, p2 and p'2, the initial and final momenta of the 

four paths are set as follows: 71 (pi < p'l), Ji{Pi * P2), 72(^2 < P2) and 

72 (P2 ^ -Pi)- 

It is convenient to compute the integral (8) in the rotated coordinates 




and 



V2 



b' 



(9) 



It will be possible to integrate over da da' in the SP approximation. Consider 
the phase function 

$(0, b, a', b') = Sj, (pi, -p[) - Sj'^ {p'i,-P2) + iP2, -p'2) - {p'2, (10) 

where the momenta must be expressed in terms of the new coordinates. It is 
easy to show that the SP condition d^/da — d^/da' = can be rewritten as 
a relation between the initial and final positions y*'-^ of the trajectories on the 
lead: 

These equations are satisfied if 



- yt,i + j/4 + y;^ = -y;, + 4 + - j/4 = o. (ii) 



71 = 72 and j[ 



l'2, 



(12) 



and the stationary point is given by a = a' = 0. Now it is necessary to employ 
the diagonal-approximation requirement $ = 0. It limits possible trajectories to 
7i = (7i)^^ =72 = (72)^^1 as expected. Since the condition $ = is fulfilled 
for arbitrary b and b' , the integral over dbdb' must be performed exactly. In the 
original coordinates, it is the integral over the manifold 



/ Pi = P2 



p'l = p'2 



(|ft|,|p',|<M- 



(13) 



To complete the program, we compute the second derivatives of the phase 
at the SP point: 



92$ 
9^ 



^2$ 



o=a'=0 



= 0, 



a=a'=0 



dada' 



= 2 



d^S^i{Pi.-p'i) 



(14) 



a=a'=0 



dpi d{-p[) 
8), accordin; 

Finally, transforming the dbdb' integral to the dpi dp'i integral, we arrive at (5). 



The SP integration leads to the cancellation of l^^iP in (8), according to (4). 
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It is instructive to consider an analogy with the semiclassical treatment of 
normal reflection. Suppose, the scattering matrix {S2) describes the prop- 
agation of the particle before (after) the reflection and S = S2S1 is the total 
scattering matrix. Scmiclassically, 5*1 and S2 arc given as sums over classi- 
cal trajectories. If S is calculated in the SP approximation, it will be a sum 
over combined classical trajectories which include the reflection. In the case 
of Andreev billiard, the consequence of the SP approximation is the equality 
of trajectories (12), while the Andreev reflection is required by the diagonal 
approximation. 

2.2.2 Quantum corrections 

There are two types of quantum corrections to (5): diagonal and off-diagonal. 
(The situation here is analogous to the calculation of the quantum corrections to 
reflection in a normal billiard [9].) The off-diagonal corrections result from the 
Sieber-Richter (SR) [9, 13] pairs of trajectories. Such trajectories are exponen- 
tially close in the phase space (up to a time reversal), apart from a small region 
where one of them has a self-crossing and the other has an anticrossing. Their 
actions differ by a small amount A<S(£:) depending on the crossing angle e. With- 
out the diagonality requirement $ = 0, a more general condition (12) should 
be used. If 71 = 72 has a crossing (anticrossing), it can be paired with j[ = 72 
which begins and ends exponentially close to (71)"^, but has an anticrossing 
(crossing). Their phase difference $ = ±lAS{e-y^) (see appendix A.l for / > 2) 
is small compared to their actions. Hence, the momentum integration, as above, 
yields the SR correction to (5) 

{Tr[S{E)S*i-E)f)^^ = 

2 ^ ^ \Ay\^ eyi^ l^2lErA cos 

n,m=l {7(n,m)} 

where the index {7} runs over all self- crossings of the path 7. Performing the 
summation over the crossings by the standard procedure [9, 13-15], we arrive at 
the final result 

{T,[S{E)S*{-E)]')^^ = -j^ ^ ^ \A,\'T,e^^(^2lETA (16) 

'^^'^ n,m=l 7(n,m) ^ ^ 

depending on the average escape time T^sc = ttA/wvf, where vp is the Fermi 
velocity (appendix B). 

The diagonal quantum correction is 

(Tr [S{E)S*{-E)]%^^ = T E E 1-^71' exp U^IET,) . (17) 

n=l ■y{n,n) 

For I = 1, it is readily derived using equation (2). Namely, in (Tr S'S'*) = 
Tlinm 7' (■ ■ ■) terms with n = m and 7 = 7' are considered. Since the ac- 
tions of the identical paths cancel in the phase, these terms survive the averaging 
over kp and yield the above result. Note that the terms with the time-reversed 
paths 7 = 7'"^ enter the classical part (5). In the case / > 1, the preceding 



(15) 
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integration procedure needs to be modified. This is done in appendix A. 2. It 
is worth mentioning that both the SR and diagonal corrections are of the next 
order in to the classical contribution. 

For the following calculation, it will be helpful to have the sums over the 
trajectories in (16) and (17) transformed into integrals over time. This can be 
achieved by using the sum rule for chaotic billiards [9]: 



(18) 



7(n,m) 



where P{T) = (T^sc) ^ exp(— T/7^sc) is the survival probability. Adding the SR 
and diagonal contributions together, we obtain the total quantum correction 



{Tv[SiE)S*i-E)]^) 



1 



quant 



dTP{T)(l-^]oxp{^2lET]. (19) 



2.3 Quantum correction to the density of states 

The quantum correction to the DoS, originating from the SR and diagonal 
contributions, 



ET 



cos ■ 



(20) 



is found by substituting the corrections to the traces (19) in equation (1). In 
deriving (20), the sum over I was computed as follows: 



ReE 



1=1 



^ exp ( ^^lET^ = -Re In 



In 2 



1 + exp ( ^2ET] 



cos ■ 



ET 



(21) 



Important conclusions can be drawn from (20) already in the limit 
E/ETh 1, where Exh = ^/27^sc is the Thouless energy. In this case, one 
finds 

ln2- 



5diE) 



TrErh 



3 E Y 
2 Et^J 



We recall that the classical part of the DoS [4, 16] 

coshx _ TrE'Th 



dci{E) = dox^ 



sinh^ X ' 



E 



(22) 



(23) 



becomes exponentially small in this limit. Hence, equation (22) implies that the 
total DoS is negative at small energies. This unphysical result has two possible 
explanations: either there are other sources of quantum corrections not taken 
into account in the present work or the semiclassical theory becomes inapplicable 
near Ep, thereby reflecting the existence of the gap. Assuming the latter, we can 
roughly estimate the gap width Eg by setting dc\{Eg) + Sd{Eg) = 0. This energy 
is weakly dependent on the channel number N and reads Eg/ETh ~ 0.45, 0.34 
for N = 10, 100, respectively. These numbers are comparable with the results 
of the random-matrix theory [3] Eg/ETh ~ 0.6, as well as the full quantum 
calculations [4]. 
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3 Conductance 



In this section, we derive semiclassical formulae for the conductance G of a 
chaotic billiard having one normal (N) and one superconducting (S) lead. It was 
shown [6, 7] that G = (4e^//i) Ra, where Ra is the Andreev-reflection coefficient 
defined in the introduction. It can be expressed [10] in terms of the electron 
transmission matrices fgN and between the leads as 

Ra = T^ [tsN^NS (2 - isNiNs)"'] ^ (24) 

(it is assumed that E = Q). Expanding the denominators in the Taylor series, 
we obtain 

oo 

Ra=Y. 2-('+'')'IV(tsNt;,s)'+''. (25) 
i,i'=i 

which is a convenient starting point for the semiclassical treatment. 

The semiclassical expressions for tsN and iNS are of the form (2) where the 
trajectories 7 now connect the respective leads. The calculation of the traces 
in equation (25) (averaged over kp) is completely analogous to that of the 
preceding section. For the classical contribution we find 

^A,ci = EE E 1^7!' = ^01, (26) 

n=l m=l 7(n,m) 

where A?s^n is the number of channels in the leads and T^x is the classical trans- 
mission coefficient for electrons. It was taken into account that the classical 
traces are independent oi I + V [E = 0) and that EO'=i 2"^'"'"''-* = 1- This 
result supports the classical argument, according to which all trajectories that 
reach the superconducting lead will contribute to Ra,c\- The quantum correction 
due to the SR pairs (there is no diagonal correction) becomes 

Ra,sk = TsR ^ [(/ + r)2'+''] = TsR In 2, (27) 
i,i'=i 

where Tsr, < is the standard SR correction to the electron transmission co- 
efficient [9]. Thus, the weak- localisation correction is smaller in an Andreev 
billiard than in the normal billiard of the same shape. 



4 Summary and conclusions 

We presented a new method which allows us to calculate traces of semiclassi- 
cal scattering matrices in the Andreev billiards. The method was applied to 
determine the density of states and conductance in the chaotic cavities. The 
classical contribution to these quantities was found to be in agreement with 
the existing theories. The current framework made it possible to compute the 
quantum corrections to the classical results. We have shown that the weak- 
localisation correction in the two-probe geometry is reduced, compared to the 
normal billiard of the same shape. In the closed cavity, the quantum corrections 
make the density of states negative within a small distance (of the order of the 



7 



Thouless energy) from the Fermi level. If this property is a signature of the 
gap in the density of states discovered by other methods, it will be important 
to understand the failure of the semiclassical theory close to the Fermi energy. 
Alternatively, other, so far unknown, types of quantum corrections could com- 
pensate the negative value. 
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A Derivations for arbitrary / 

A.l Classical contribution 

We start with the generalization of integral (8) 



Tr [S{E)S* (-£)]' = ( — ) / dpi • • • dpidp'i 

X A^,A;, ■ ■ ■ A^,A;, exp U$ j exp U^lETj , (28) 



71 .7i. 
lull 



where the paths, in terms of their initial and final momenta, are ^j{pj < p'j) 

and 7^ (p^ < pj+i) {j = 1, . . . J; I + 1 = 1), ^ = Y!']^-^ {s^^ ~ ^i'^ 

T = {2l)~^ X^j=i {j^ij + '^I'j)' ^® transform to the new variables 



/ ai 










/ a[ 


\ 










and 






ai-i 














V b 


) 








J 



in such a way that b and 6' are within the manifold 

{\Pi\,\p'i\ < kp), 



Pi 



Pi 



(29) 



Pi 
P'l- 



--Pi=P 
P'l^f 



while ai and are transverse to it. This means that C is an orthogonal matrix 

which has the property ^ij = (i = 1, 

choose this matrix in the form 



(30) 
itrix 

,1—1). It is convenient to 





\ 7b(i. 


-1, 


0, 


0, 


0, .. 


., 0, 


0) 






1, 


-2, 


0, 


0, .. 


., 0, 


0) 


c = 




1, 


1, 


-3, 


0, .. 


., 0, 


0) 






1, 


1, 


1, 


1, .. 


1, 


-{I -I)) 






1, 


1, 


1, 


1, .. 


• , 1, 


1) 



(31) 
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where the rows are niultiphcd by the factors in front of them. 
The derivatives of the phase function 



dak 



(32) 



and similar for 9$ / da'f^ . vanish if 71 = • • • = 7; = 7 and 7^ = • • • = 7; = 7' . In 
the diagonal approximation, we have $ = and 7 = (7')"^- In the SR quantum 
correction, 7 and (7')"^ form the SR pair, and the phase is $ = ±ZA<S(£y). 

Wc proceed in the diagonal approximation. The non-vanishing second deri- 
vatives of $ at the SP point are 



92$ 



daida'i. 



=a' =0 



dpd{—jf) 



■D. 



(33) 



where wc introduced the {I — 1) X {I — 1) matrix Dik = Sik — C^i^j+iCj^. 

Equation (31) yields an explicit expression Di^ = Dik/\/i{i + l)k{k + 1) in 
terms of the matrix 



D : 



1-2 + 1 
1 
1 



-1 • 3 
2-3 + 1 
1 





-2 - 4 
3-4 + 1 






-3 - 5 



1 1 1 1 1 ■■■ 1 (l-2)(/-l) + l -{1-2)1 

\ I I I I I ■ I I J 

(34) 

Its determinant det D = (Z!)^ can be computed by adding the first column to the 
second column, then adding the resulting second column to the third column 
and so on. Hence, the Hessian of $ is given by 



det 



dai da'. 



(j)l±_\ ( 

y da'^dak J \ 

To complete the SP integration 



da'- da'. 



d^Sj{p, -p' 
dpd{—p') 



2(1-1) 



I'. 



(35) 



/ 



I^^^^exp -cl> + z-sgncl> )(- 



(36) 

we find the difference of the number of positive and negative eigenvalues of 

A 

sgn$" = 0. This can be shown as follows: if an eigenvector of ' 



with an eigenvalue A is 



^ ^ , then ^ ^ is an eigenvector with the eigen- 
value —A (here A and B are the square matrices and a and /? are the columns 
of equal size) . After the last two integrations over db db' along the manifold (30) 
we end up with equation (5). 



A. 2 Diagonal quantum corrections 

Let us repeat the SP integration over 0^=1 dajda'^ in appendix A.l with the 
exception that now 7 = 7', instead of 7 = (7')"^) is taken at the stationary 
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point. This choice is compatible with the diagonahty condition $ = and is 
possible on the manifold p = p' . It is convenient to perform the remaining 
integrations over dbdb' in the transformed coordinates b± = l^^^^{b± 6')/\/2 
such that $ = is fulfilled on the line 6_ = 0. The subsequent averaging over 
kp is expected to pick up the contribution of the neighbourhood of this line 
(diagonal approximation). Therefore, we can expand $ « {d^/db-)b-, where 
d^/db- = {yij - yt^) / ^/2. This leads to the result 



/ 




rf6_exp( ) ). (37) 



Note that this formula includes two classes of trajectories, 7(p <~- —p) and 
7(p^p), and, correspondingly, two signs of j/^ in the second exponent. The 
latter class is taken into account if one starts from the equation (28) with 

iAP] ^P'j) and iAp'. ^ Pj 



-1) 



One can avoid explicit calculation of equation (37) by comparing it with the 
result (17) for / = 1, which was derived independently. It is then quite obvious 
that (17) is valid for arbitrary I. 



B Summation over the Sieber-Richter pairs 



The purpose of this section is to fill in the steps between equations (15) and (16). 
For a trajectory 7, the average number of self-crossings with the crossing angle 
between s and £ + ck is [9] 



Px{e;T,)de: 



ttA 



sme 



1 _ 2,rmi„(£) 



de, 



(38) 



where 7^in(£) is the size of the crossing region (logarithmically dependent on e). 
We took into account the correction put forward in [15], according to which the 
%nin{£) contribution in (38) is twice as large as was previously suggested in [9]. 
The terms of higher order in 7^in(e)/X^ are neglected. It is important to keep 
the contribution proportional to 7min(e)i since the subsequent integration over 
e would eliminate the leading term [13]. 

With the help of (38), the sum over the self-crossings in (15) can be reduced 
to the sum over the paths as ^[^(^n,m)} (' ' ' ) ^ / J2-y{n,m;£) 
where the index 7(n, m; e) runs over all paths having a self-crossing of angle e. 
Let us, for a moment, transform the latter sum to an integral using the sum 
rule (18). It was noticed by the authors of [15] that the sum rule has a correction 
linear in 7^in(e)- It can be obtained by shifting T i— » T — 7^in(£) in the right- 
hand side of (18). Thus, there are two contributions of the first order in 7^in(£): 
the one coming from the leading-order term in (38) and the first-order correction 
to the sum rule and the other resulting from the '7min(£) term in (38) and the 
uncorrected sum rule. Integrating over T explicitly, one can show that the 
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former contribution is two times smaller than the latter and has the opposite 
sign. 

The preceding argument allows, in effect, us to consider the half of the first- 
order term in (38) and transform X/{7(n,m)} (' ' ' ) ^ X^7(n,m) 
The e integral was computed in [13] and reads 



4 r 



ds cos 



IAS (c) 



sin(£)T^i„(£)«^^^^. (39) 



h 

Then equation (16) would follow. 
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